/**
 * @file   InterpolationSolver.h
 * @author WinterMelonCDY <wintermeloncdy@wintermeloncdy-VirtualBox>
 * @date   Tue Oct 25 10:09:47 2022
 * 
 * @brief  
 * 
 * 
 */

#include <iostream>
#include <cmath>
#include <algorithm>
#include <limits>
#include <vector>

#define precision 1e-7 
//[题目A要求]Implement the Newton formula in a subroutine that produces the value of the interpolation polynomial pn(f; x0, x1, . . . , xn; x) at any real x, where n ∈ N+, xi’s are distinct, and f is a function assumed to be available in the form of a subroutine.

//定义函数类
class Function
{
public:
    virtual double operator()(double _x) = 0;
};

//定义插值类
class Interpolation
{
public:
    virtual double solve(double _x) = 0; 
};

//Newton插值多项式继承函数类
class Newton_polynomial : public Function
{
private:
    std::vector<double> C, P;
    int n;
public:
    Newton_polynomial(std::vector<double> _C, std::vector<double> _P): C(_C), P(_P) {}
    double operator()(double _x)
    {
        double y = C[0];
        n = P.size();
        for (int i = 1; i < n; i++)
        {
            double delta_y = C[i];
            for (int j = 0; j < i; j++)
            {
                delta_y = delta_y*(_x - P[j]);
            }
            y = y + delta_y;
        }
        return y;
    }
};

//Hermite插值多项式继承函数类
class Hermite_polynomial : public Function
{
private:
    std::vector<double> C, L;
    int n;
public:
    Hermite_polynomial(std::vector<double> _C, std::vector<double> _L): C(_C), L(_L) {}
    double operator()(double _x)
    {
        double y = C[0];
        n = L.size();
        for (int i = 1; i < n; i++)
        {
            double delta_y = C[i];
            for (int j = 0; j < i; j++)
            {
                delta_y = delta_y*(_x - L[j]);
            }
            y = y + delta_y;
        }
        return y;
    }
};

//Newton插值公式的运用，最后利用差商求出Newton插值多项式系数
class Newton_Formula : public Interpolation
{
private:
    std::vector<double> X, Coef;
    Function &f;
    int n;
public:
    Newton_Formula(Function &_f, std::vector<double> _X): f(_f), X(_X) 
    {
        n = X.size();
        double table[n][n];
        for (int i = 0; i < n; i++)
        {
            table[i][0] = f(X[i]);
        }
        for (int i = 1; i < n; i++)
        {
            for (int j = i; j < n; j++)
            {
                table[j][i] = (table[j][i-1] - table[j-1][i-1])/(X[j] - X[j-i]);
            }
        }
        for (int i = 0; i < n; i++)
        {
            Coef.push_back(table[i][i]);
        }
    }
    
    //solve函数返回Newton插值多项式任意一点对应的多项式值
    double solve(double _x)
    {
        Newton_polynomial Pn(Coef, X);
        return Pn(_x);
    }

    //Get_coef函数返回Newton插值多项式的相应系数
    std::vector<double> Get_coef()
    {
        return Coef;
    }
};

//Hermite插值公式的运用，最后利用差商求出Hermite插值多项式系数
class Hermite_Formula : public Interpolation
{
private:
    std::vector<double> X, Label, Coef;
    int n;
public:
    Hermite_Formula(std::vector<double> _Label, std::vector<double> _X): Label(_Label), X(_X)
    {
        n = X.size();
        double table[n][n] = {{0}};
        int count = 1;
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < count; j++)
            {
                table[i][j] = X[i+j+1-count];
            }
            if (Label[i] == Label[i+1])
                count = count + 1;
            else
                count = 1;
        }
        for (int i = 1; i < n; i++)
        {
            for (int j = i; j < n; j++)
            {
                if (table[j][i] == 0)
                    table[j][i] = (table[j][i-1] - table[j-1][i-1])/(Label[j] - Label[j-i]);
            }
        }
        for (int i = 0; i < n; i++)
        {
            Coef.push_back(table[i][i]);
        }
    }
    
    //solve函数返回Hermite插值多项式任意一点对应的多项式值
    double solve(double _x)
    {
        Newton_polynomial Pn(Coef, Label);
        return Pn(_x);
    }
    
    //solvediff函数返回插值多项式在_x点处的导数
    double solvediff(double _x)
    {
        Newton_polynomial Pn(Coef, Label);
        return (Pn(_x + precision) - Pn(_x - precision))/(2*precision);
    }

    //Get_coef函数返回Hermite插值多项式的相应系数
    std::vector<double> Get_coef()
    {
        return Coef;
    }
};
